Self-interaction correction with Wannier functions 
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We describe the behavior of the Perdew-Zunger self-interaction-corrected local density approxi- 
mation (SIC-LDA) functional when implemented in a plane-wave pseudopotential formalism with 
Wannier functions. Prototypical semiconductors and wide-bandgap oxides show a large overcorrec- 
tion of the LDA bandgap. Application to transition-metal oxides and elements with d electrons is 
hindered by a serious breaking of the spherical symmetry, which appears even in a closed shell free 
atom. Our results indicate that, when all spherical approximations are lifted, the general applicabil- 
ity of orbital-dependent potentials is very limited and should be reconsidered in favor of rotationally 
invariant functionals. 
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I. INTRODUCTION 

Within the Kohn-Sham approach to density functional 
theory, the total energy of a many-electron system of 
density p = p-f + pi is generally decomposed into four 
terms: 



E KS = T s + E H [p] + E ext + E xc [ PhPl }. 



(1) 



These terms describe the kinetic energy of the fictitious 
set of one-particle orbitals T s , the Hartree energy Eh, 
the energy due to the interaction with the external po- 
tential E ext , and the unknown exchange and correlation 
energy E xc ; the latter is commonly approximated by lo- 
cal or semi-local functionals of the density such as the 
local density (LDA) or generalized gradient (GGA) ap- 
proximations. 

Despite the innumerable successes of the LDA and 
GGA, some serious drawbacks exist that prevent the ap- 
plicability of these methods to a wider range of materials 
and phenomena. Situations in which these standard func- 
tionals lead to qualitatively incorrect physics include the 
erroneous prediction of metallicity for magnetic transi- 
tion metal oxides, an inability to localize defect states in 
solids 1 and unpaired electrons in water 2 , qualitatively in- 
correct metallic transport for single- molecule junctions 3 , 
inaccurate redox potentials and charge-transfer reac- 
tions^, and unphysical fractionally charged fragments in 
the molecular dissociation limit 5 . These failures can be 
traced, at least in part, to the self-interaction error (SIE), 
which is the spurious interaction of an electron with its 
own Hartree and exchange-correlation potential. Indeed, 
in the case of one-electron systems such as the ground 
state of the hydrogen atom, Eh and E xc should cancel 
exactly: 

Eh[Pis] + E xc [p ls ,0] should equal 0, 

but this condition is not fulfilled by approximate 
exchange-correlation functionals such as LDA or GGA. 
While in a many-electron system the notion of self- 
interaction is less clear cut, it is commonly accepted that 



this same mechanism affects the behavior of strongly lo- 
calized, atomic-like orbitals, such as d states in tran- 
sition metal compounds, by suppressing or mistreating 
on-site Coulomb interactions. The considerable funda- 
mental and technological interest in d-electron systems 
such as high-T c superconductors and colossal magnetore- 
sistive manganites provides a compelling incentive for im- 
plementing appropriate SIE- free density functional meth- 
ods. Interestingly, the SIE has a relatively minor impact 
on total energies, but it strongly affects the eigenvalues 
of the Kohn-Sham Hamiltonian. In particular, the en- 
ergy eigenvalue associated with the highest occupied or- 
bital usually shows a strong departure from the ionization 
potential, while it should match it exactly within exact 
DFT 6 . 

Attempts to correct the self interaction error can be 
traced back to the seminal paper by Perdew and Zunger 
(PZ) 7 , who defined the self- interaction corrected (SIC) 
exchange-correlation energy, E XI ! C , as 

E S x l C = ^r° x [PT,P|]-E PhIPo*] + E*™™[p aa , 0]) . 



(2) 

Here E^ pTOX [p^, pjj is the approximate (for example 
LDA or GGA) exchange-correlation energy, and the term 
within the summation is the self-interaction energy of an 
electron in orbital a with spin a; Enlpaa-] is the self- 
Coulomb part and E x Pprox [p at7 , 0] is the self exchange- 
correlation part. For isolated atoms, this approach 
yielded Hamiltonian eigenvalues which were in surpris- 
ingly good agreement with experimental removal ener- 
gies. 

These successes motivated a considerable subsequent 
effort to incorporate PZ self-interaction corrections in cal- 
culations for solids. Unfortunately, however, direct im- 
plementation of the PZ functional in extended systems 
has proved to be technically non-trivial. The main issue 
arises from the fact that the SIC-LDA functional, un- 
like standard Kohn-Sham theories, is not invariant with 
respect to a unitary transformation of the occupied man- 
ifold; in particular, the SIC vanishes for extended Bloch 
wavefunctions. Therefore, the first challenge of any im- 
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plementation is to devise a general and physically sound 
criterion for the choice of this unitary transformation, 
which yields a set of "local orbitals" (LO), as opposed to 
the "canonical orbitals" (CO) which are the usual eigen- 
states of the Hamiltonian with Bloch periodicity 

Soon after the initial work by Perdew and Zunger, 
Heaton Harrison and Lin (HHL) recognized that Wan- 
nier functions provide an ideal basis for describing the 
localized-delocalized duality of electrons in the full- 
SIC Hamiltonian^; by implementing SIC-LDA within a 
LCAO basis set HHL found considerable improvement in 
the solid Ar and LiCl bandstructures. An appealing as- 
pect of the HHL approach is the introduction of a unified 
Hamiltonian by means of band projections. This strategy 
removes the orbital dependence of the SIC Hamiltonian 
and allows for the calculation of all SIC-LDA eigenvalues 
for a given k point by one single matrix diagonalization. 
Furthermore, HHL defined the unitary transformation 
between the Wannier and Bloch representation as the 
one yielding the variational minimum of the SIC-LDA 
functional within the usual orbital orthonormality con- 
straints. Pederson, Heaton and Lin 9 later demonstrated 
that the so-called "localization condition" is then fulfilled 
by the localized orbitals 4> a and their associated SIC po- 
tentials 5V a : 

(4> a \SV a -6V^ ). (3) 

This means that the Lagrange multiplier matrix enforc- 
ing the orthonormality condition is Hermitian, and it can 
be indeed diagonalized to obtain the SIC eigenvalue spec- 
trum, together with a set of eigenvectors that correspond 
to the Bloch-periodic COs. 

Svane and GunnarssoniS (SG) and later Szotek, Tem- 
merman and Winter— (STW) applied a fully self- 
consistent SIC-LDA method to extended systems within 
an LMTO-ASA (linear muffin tin orbital - atomic sphere 
approximation) implementation, obtaining remarkable 
results for both d- and /-electron materials. A major 
pitfall of the SIC functional is that it allows for multi- 
ple local minima, one of these being the trivial solution 
where all the LOs are Bloch- like (i.e. no SIC is applied), 
and another obvious one being the one where all the LOs 
are Wannier-like; intermediate (mixed) choices also ex- 
ist, where some of the LOs are Wannier-like, and others 
keep their itinerant character. SG and STW proposed 
choosing the solution with the lowest total energy (which 
corresponds to the absolute minimum of the SIC func- 
tional, and is consistent with the variational character of 
the localization procedure). Based on this choice, phase 
transitions are sometimes observed as a function of exter- 
nal parameters (e.g. volume) in which the SIC contribu- 
tion for a given band becomes positive (or negative); this 
crossover between SIC and no SIC is rationalized in terms 
of a physically appealing realization of Mott transitions 
(which are driven by the competition between bandwidth 
and on-site Coulomb repulsion). 

Some fundamental problems with the use of a partially 
Bloch- like and localized solution have been pointed out, 



however, by Arai and Fujiwara 12 (AF). First, the pres- 
ence of "delocalized" bands to which no SIC is applied 
leads to a severe size-consistency problem when the ex- 
tended solid is considered as the thermodynamic limit 
of an increasingly large cluster where SIC is unavoid- 
ably finite^. Even in regions where the SIC energies 
are slightly positive the SIC potentials remain strongly 
attractive, so when the cluster volume is increased to 
the thermodynamic limit strong and unphysical changes 
in the eigenvalue spectrum must be expected. Second, 
the sign of the SIC energy (and hence whether or not 
the orbitals are treated as localized) is sensitive to de- 
tails such as the parameterization of the LSDA. Since 
the main aim of the SIC method is to correct for the self- 
exchange error, qualitative differences in the electronic 
ground state determined solely by minor details of the 
correlation functional are, again, physically hard to jus- 
tify. 

Interestingly, AF also discussed the consequences of 
the sphericalization of the SIC potential, which is rou- 
tinely performed (see SG and STW) within LMTO im- 
plementations and was also adopted in the early works of 
HHL 8 . While a significant impact on the bandstructure 
of solids and unphysical energy splittings within other- 
wise crystal-symmetric multiplets were found, AF con- 
cluded that the overall consequences of this approxima- 
tion were relatively unimportant. 

It has been shown recently for a wide range of atoms 
that orbitals with different angular momenta are allowed 
to mix upon lifting the spherical approximation^ 3 -. For 
example 3s and 3p states in Ar mix to yield four tetra- 
hedrally symmetrical sp 3 hybrids, analogous to the max- 
imally localized Wannier functions in sp compounds^. 
As a consequence of this mixing, the SIC energy be- 
comes negative for all bands, while it is generally pos- 
itive for pure p states; furthermore, eigenvalue shifts of 
the order of 1 eV occur, and in general the agreement 
with experiments tends to worsen 13 . These atomic re- 
sults suggest that, while the sphericalization of the SIC 
potential itself has a very minor direct impact, in agree- 
ment with the conclusions of Refs. H and [IH, it might 
well have a dramatic indirect impact, by preventing the 
true variational minimum of the SIC functional from 
being found. In particular, if the total energy crite- 
rion for the selection of the localized/delocalized bands 
is enforced, the artificial suppression of interband mix- 
ing might lead to erroneous choices, and qualitatively 
incorrect electronic ground states; for example, oxygen 
2p states, that are considered to be itinerant within the 
spherical approximatio n 10 ' 11 can become localized once 
mixing with oxygen 2s states is allowed. 

It is apparent from the above analysis that two main 
issues affect self-interaction corrected calculations for ex- 
tended solids, i) the existence of multiple local min- 
ima and ii) the validity of the spherical approxima- 
tion. In this work we address both issues by using mod- 
ern Wannier-function theory^ and a plane-wave norm- 
conserving pseudopotential formalism. By testing our 
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method on simple atomic systems we first demonstrate 
that, if used with care, the pseudopotential approxima- 
tion introduces a negligible error with respect to the 
most accurate all-electron results available to date 13 ; 
this provides a strong validation of our results and con- 
tributes to putting the full-SIC formalism onto a solid 
implementation-independent technical footing. In par- 
ticular, in agreement with Ref. Il3l . we find that the 
spherical approximation has an important impact on the 
eigenvalue spectrum of solids, often significantly worsen- 
ing the agreement with experimental spectroscopic data. 
We further find that, within our spherically- unrestricted 
SIC-LDA, the fully localized solution is always the varia- 
tional electronic ground state, even in bulk Si where the 
valence electrons are usually considered as being itin- 
erant; this result suggests that some caution must be 
taken when interpretating the localized/delocalized SIC 
crossover in terms of a Mott transition, since it might 
be an artifact of the numerical approximations used. Fi- 
nally, in our implementation two further pitfalls of the 
SIC-LDA method emerge, which were so far overlooked 
in the literature: i) the giant overcorrection of the elec- 
tronic band gaps in solids and ii) the unphysical break- 
down of crystal point symmetry, which is especially seri- 
ous in ci-electron systems. We rationalize these effects in 
terms of, respectively, lack of proper treatment of dielec- 
tric screening and the rotational non-invariance of the 
method. Our results provide useful guidelines for further 
research in the quest for an improved density functional, 
and also a benchmark against which approximate flavors 
of SIG±&±£i±£ can be tested and validated. 

The remainder of this work is organized as follows. In 
Sec. |TT]we give a detailed overview of the SIC technique 
we use in this work. In Sec. IIIII we present our results: 
First we validate our method by performing some tests 
on simple atomic systems, then we apply SIC to simple, 
prototypical solids (Ar, Si, MgO) and finally, we analyze 
the performance of SIC for d-electron systems. In Sec. lIVI 
we discuss these results in light of previously reported 
studies and analyze their impact for future methodologi- 
cal development. In Sec. [V] we summarize and conclude. 
The Appendix presents an analysis of Boys orbitals in 
d-electron spherically symmetric atoms; this analysis ex- 
tends the work by Posternakii for sp elements, and shows 
that orbital-dependent functionals tend to unphysically 
break the symmetry of closed-shell atoms. 



alized Bloch orbitals ^nk(r) (which are not necessarily 
eigenstates of the Hamiltonian) can be written in terms 
of the cell-periodic functions u„k(r): 

^nk(r) = e lk - r u„ k (r) ; 

the latter can be represented in reciprocal space as fol- 
lows: 

Unk(G) = —7= ( e^ Gr u„ k (r)dr . 

V" J cell 

The reciprocal lattice vectors b of the Born- von Karman 
supercell can be written in terms of G and the fc-point 
mesh: 

b = G + k , 

and the Wannier function associated with the band n and 
the lattice site R in reciprocal space is: 



N 



(4) 



We will set R = in the remainder of the paper, and thus 
focus on the minimal set of Wannier functions which is 
necessary to describe the solid; also we will introduce a 
spin index a. Upon Fourier transformation one obtains 
the Wannier functions w na (r) in real space and their as- 
sociated charge density distributions /w(r) = |u>„ CT (r)| 2 . 
The self-interaction energy of the system, which needs 
to be added to the LDA (or GGA) total energy, is then 
given by Eqn. O 

In order to minimize the SIC-LDA functional we need 
to calculate gradients of the SIC energy with respect to 
the wavefunction plane-wave coefficients. We start by 
calculating the gradients of the SIC energy with respect 
to the Wannier functions, which can be written as: 



SE 



sic 



= Vl IC (r)w na (r) 



Here V^J C is the SIC (Hartree plus exchange and correla- 
tion) potential generated by the Wannier density p nry (r) . 
The state-dependent potential V^J can be recast into a 
unified operator by using band projections: 



V 



sic 



(5) 



II. METHOD 

Before presenting our method we briefly summarize 
some basic notations and conventions that will be use- 
ful in the derivation (for a more extensive treatment see 
Ref. [l^ and references therein). We assume a Born- von 
Karman supercell of volume £IbvK = NQ, where N is 
the total number of k points arranged on a regular three- 
dimensional mesh and £1 is the volume of the primitive 
cell used to represent the periodic crystal. The gener- 



In general, V SIC has nonzero Hermitian and anti- 
Hermitian components: 



V 



SIC-H 



V 



SIC-A 



(V 
(V 



sic 



sic 



ySIC^ 



(G) 

(7) 



When applied to the electronic wavefunctions, the anti- 
Hermitian part y SIC ~ A produces a unitary mixing 
within the occupied manifold and is the signature of the 
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rotational non-invariance of the SIC-LDA functional; no 
such term exists in standard Kohn-Sham theories. The 
Hermitian term, on the other hand, evolves the electronic 
subsystem in a direction which is perpendicular to the 
occupied subspace, and can be treated as an additional 
term to be added to the Kohn-Sham Hamiltonian. For 
reasons of transparency and in order to have better con- 
trol over the minimization process we decided to separate 
the two tasks into two nested loops. 

In an inner loop, we constrain the update of the wave- 
functions to a unitary mixing within the occupied mani- 
fold: 



and we seek the set of unitary matrices CA k ) that yields 
the minimum value of the SIC energy (the standard 
Kohn-Sham energy is invariant with respect to this trans- 
formation) . This operation is formally analogous^ to the 
calculation of the maximally localized Wannier functions 
for a set of entangled energy bands^, except that, in- 
stead of minimizing the quadratic spread, here we need 
to enforce the representation with the minimum value of 
the Perdew-Zunger self-interaction. In particular, the ro- 
tation matrices are obtained by adding an infinitesimal 
anti-Hcrmitian matrix to the identity: 



1 + dW {k) 



the variation of the SIC-LDA functional with respect to 
this transformation is provided by y SIC '- A : 



dE 



sic 



SIC-A 



\lpnk 



It is then easy to show that the stationarity of the func- 
tional (zero gradient) implies: 

(w m \V„ - V m \w n ) = , 

which is the usual "localization condition"—. 

In the outer loop we add to the LDA Hamiltonian the 
Hermitian part of the SIC operator: 



H 



SIC-LDA 



H 



LDA 



V 



SIC-H 



which is now identical to the full SIC operator since 
the anti-Hermitian part vanishes within the subspace 
spanned by the occupied bands. We then take stan- 
dard electronic steps until the ground state is reached. 
At self-consistency, the eigenvalues of this SIC Hamilto- 
nian formally agree with the eigenvalues of the Lagrange 
multiplier matrix that can be obtained within direct min- 
imization techniques 2 ^. 

We note that particular care must be taken in the cor- 
rect evaluation of the sclf-Hartree energy and potential 
of the Wannier charges, since the periodic boundary con- 
ditions induce some spurious long-range interactions be- 
tween the localized charge distributions. While some au- 
thors 22 have proposed to truncate the 1/r Coulomb inter- 
action to eliminate the divergence for G = 0, in this work 



we use the standard approach of introducing a uniform 
background charge to neutralize the system: 



e?\p] = J 



4tt 



2ft 



BvK 



E 



\pM s 
b 2 



(8) 



For a cubic BvK cell of dimension L, the error due to 
the use of periodic boundary conditions can be corrected 
up to the order 0(L~ 5 ) by the following term 2 ^: 



^=^ + |a£/ 3 ^(r-ro)r 2 



(9) 



where ro is the center of the Wannier function charge. We 
use the same form for the Hartree potential, which is the 
analytic derivative of the above term with respect to the 
charge density p n ; the relationship E H — 1/2 fV H pd 3 r 
is exactly respected. 

For both the inner and the outer loops we use a 
damped-dynamics minimization algorithm. For the for- 
mer, we checked the internal consistency of the imple- 
mentation by taking a frictionless run; the mathematical 
constant of motion was conserved within machine preci- 
sion. The latter procedure is the standard Car-Parrinello 
approach with the addition of the Hermitian SIC oper- 
ator. The method was implemented in an "in-house" 
electronic structure code. For all our tests we used a 
cubic BvK supercell, the local density approximation 
and norm-conserving pseudopotentials 24 . The atomic 
tests were performed by using a T-point only sampling 
of the Brillouin zone and a large supercell; the above al- 
gorithm did not require any modifications, since it was 
constructed to be invariant with respect to Brillouin-zone 
folding, and hence the T-point only calculations are just 
a special case. This flexibility provides an appealing link 
between isolated atoms and solids, which can be treated 
on the same footing with the exact same computational 
parameters and pseudopotentials. 



III. RESULTS 

A. Test: sp atoms, Be and Ar 

In order to check the reliability of our method we first 
apply our SIC-LDA functional to simple isolated atoms, 
with only s and p valence electrons. For consistency with 
the bulk solids, we perform these tests with the same 
plane-wave electronic structure code by placing the indi- 
vidual atoms in a large cubic cell of dimension cto = 10 A. 
For Be we obtain the eigenvalue C2 S = -9.2 eV (LDA=- 
5.6 eV), compared to -9.1 eV recently obtained with an 
all-electron SIC-LDA formalism 13 . For Ar we find the 
same level of agreement with the all-electron SIC-LDA 
result: 63 P = -16.8 eV in our calculation (LDA=-10.4 
eV), compared to -16.7 eV (all-electron). It is reassuring 
to see that the pseudopotential frozen-core approxima- 
tion, together with the adoption of periodic-boundary 
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conditions, has negligible influence on the accuracy of 
SIC-LDA in atoms, with an error of about 0.1 eV in 
a contribution that amounts to 3-6 eV. This favorable 
agreement between two very different electronic structure 
methods stems from the fact that the formalism (global 
minimum of the SIC-LDA functional with orthogonality 
constraints and no spherical averaging) is the same. 

We note that, while the Be example is trivial (there 
is only one spherically symmetric s orbital, and no opti- 
mization of the "rotational" internal degrees of freedom is 
necessary) , the Ar case has a more interesting solution in 
that the localized orbitals, just like the Boys orbitals, are 
four sp 3 hybrids with tetrahedral symmetry. The mixing 
of s and p orbitals is only allowed when the spherical ap- 
proximation is lifted, and has dramatic consequences on 
orbital eigenvalues^. Indeed, if we artificially suppress 
this mixing, we obtain e^p = -15.6 eV, a value which 
is quite close to the original Perdew-Zunger work (c3 P 
= -15.8 eV) and to the experimental ionization energy 
IE=15.8 eV. Even if the use of the spherical approxi- 
mation tends to bring atomic eigenvalues in much bet- 
ter agreement with the experimental spectroscopic data, 
this procedure is ill-defined for solids and molecules and 
therefore cannot be used as an ingredient for a general 
electronic structure method; for this reason we caution 
against its use, as did the authors of Ref. Il3l . 

We also note that, by suppressing the s — p intermix- 
ing, the SIC energy associated with each local orbital 
changes radically. The four symmetric sp 3 hybrids con- 
tribute 5 sp 3 = —0.45 eV per electron, while in the re- 
stricted solution the s orbital contributes S s = —0.52 eV 
and the p orbital S p — 0.15 eV. In addition to a much 
higher total energy, the sphericalized solution is charac- 
terized by a positive value for 5 P . Thus, the p states 
might be incorrectly discarded in the variational opti- 
mization procedure once the atom is embedded in a peri- 
odic lattice 10 , when in fact the solution with hybridized 
sp 3 states would have lower energy. 



B. sp Solids: Ar, Si and MgO 
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FIG. 1: Density of states for bulk Si (top) and MgO (bottom) 
calculated within the SIC and LDA approximations. The 
arrows and numerical values indicate the band gaps in the 
two approximations to be compared with the experimental 
values at the bottom right of each plot. 



Having assessed the reliability of our plane- wave pseu- 
dopotential implementation of the Perdew-Zunger SIC 
functional, we now move on to the more interesting case 
of solids. We choose as our examples face-centered-cubic 
(FCC) Ar, to make a direct link to the atomic tests re- 
ported in the previous subsection, a prototypical semi- 
conductor (Si) and an insulator (MgO); all these mate- 
rials show the typical LDA underestimation of the band 
gap. The motivation for investigating these apparently 
"simple" compounds is to better understand the behav- 
ior of the SIC functional in well-known test-case systems, 
before moving to more complex solids where the descrip- 
tion at the LDA level is highly problematic. We use a 
FCC primitive cell and a simple cubic Born- von Karman 
supercell corresponding to a5x5x5 (3x3x3) /c-point 
meshes in the cases of Si (MgO and Ar). The experimen- 



tal lattice constants are used in the Si and MgO cases 
(10.2 and 7.96 bohr respectively), while the lattice con- 
stant is progressively varied in the case of Ar from the 
experimental value to an artificially compressed state. 
Plane-wave cut-offs of 50 Ry, 20 Ry, 70 Ry are used in 
the respective cases of Ar, Si and MgO. 

The main goal of the calculations for solid Ar is to 
quantitatively evaluate the effect of SIC in the transition 
from the atomic problem to the bulk solid. In particular, 
we test the statement that "the SIC-LSD approximation 
provides a mechanism which allows the wave functions 
to localize for systems where the hopping integrals are 
small relative to the Coulomb interactions"—, by tuning 
the magnitude of the bandwidth, W . We use hydrostatic 
pressure to vary the bandwidth of the 3p band of FCC 
Ar from the experimental equilibrium volume (arj =9.9 
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TABLE I: Bandwidths, W, and band gaps, E g , calculated 
within the LDA and SIC approximations for solid Ar over a 
range of lattice constants. All energies are reported in eV. 
The last three columns show the SIC energy contribution per 
electron, Esic, the SIC correction to the band gap A 9 and 
the average of the SIC potential, {Vsic)- The lowest row lists 
the corresponding values for the isolated atom. 



a.u., Wlda=^-3 eV) to a highly compressed state (a 
=7.0 a.u., Wl,da=8.4 eV), and monitor the effect of SIC 
for these extremes together with two intermediate values; 
the results are reported in Table |TJ 

As expected, the LDA 3p bandwidth (Wlda) progres- 
sively increases as the crystal is compressed, and the elec- 
tronic gap increases as well. The same trend is respected 
in SIC-LDA, with a slightly larger bandwidth (by 0.1- 
0.5 eV) and a rather dramatic opening of the electronic 
gap with respect to the corresponding LDA results. The 
striking fact which is apparent from the Table is that 
the SIC correction to the electronic gap is practically 
independent of pressure, and amounts exactly to the cor- 
rection to the 3p orbital eigenvalue of the free Ar atom 
(6.46 eV). Even more striking is the lack of pressure de- 
pendence (within numerical error) of both the SIC energy 
contribution per electron E SIC (-0.44 to -0.46 eV), and 
the average value of the SIC potential on the correspond- 
ing Wannier function (-7.05 to -7.12 eV); this indicates 
that in this system the SIC is substantially insensitive 
to the bandwidth of the solid, and corrections are iden- 
tical to those calculated in the free atom. Most notably, 
E SIC is almost constant. Therefore, application of the 
SIC always lowers the variational energy of this system, 
and there is no crossover to a hypothetical delocalized 
solution. These results (together with the discussion of 
atomic Ar in the previous section) strongly suggest that 
the itinerant character of the oxygen 2p bands reported 
previously^ is a result of the spherical approximation 
adopted therein, rather than an intrinsic physical feature 
of the SIC-LSD method. 

Next we move to the cases of Si and MgO. In Fig. [1] we 
compare the calculated SIC and LDA densities of states 
and band gaps for both materials at the experimental 
lattice constants. In all cases the top of the valence band 
is set to eV. As in the case of solid Ar, the main effect 
of the SIC is an important stabilization of the valence 
bands compared to the unoccupied states; otherwise the 
density of states appears to be almost unaffected, apart 
from a slight increase of the bandwidth within SIC com- 
pared to the LDA ground state. The significant opening 
of the band-gap leads to a dramatic overcorrection of 



the LDA value, especially in the case of bulk Si. The 
band-gaps within SIC-LDA are respectively 4.5 eV for 
Si and 11.6 eV for MgO, compared to the LDA (exper- 
imental) values of 0.4 (1.2) eV and 4.6 (7.8) eV. This 
behavior might seem surprising at first sight, especially 
in silicon where the highly dispersive character of the va- 
lence bands leans heavily towards a delocalized (Bloch) 
description of the electrons rather than a localized one. 
Our results, however, suggest that even in Si the Wannier 
functions (which in this case are centered along the Si-Si 
covalent bonds) are localized enough to carry a signifi- 
cant SIC {Esic = -0.24 eV, (Vsic) = -4.58 eV); this 
fact further undermines the validity of the SIC-LDA as a 
theory to discriminate between band insulators and Mott 
insulators. As a further proof of the localized character 
of oxygen p bands in solids, we note that our SIC-LDA 
solution for bulk MgO shows similar behavior to that of 
the Ar crystal, in that four sp 3 hybrids are formed, each 
with decidedly negative values of Esic = —0.47 eV and 
(Vsic) = -7.80 eV. 

We are aware of three previous SIC calculations for 
these materials. HHL 8 found a correction of 6 eV for 
the bandgap of solid Ar, which is fairly close to our re- 
sult in spite of use of the atomic orbital and spherical 
approximations in the earlier work. FCC Ar was also 
investigated by Szotek, Temmerman and Winter—, who 
found an increase of the bandgap of 5.1 eV only, which is 
in better agreement with our results for the spherically 
restricted atom. Bulk Si was studied within an approx- 
imate bond self-interaction correction by Hatsugai and 
Fujiwara2£, who found a very favorable agreement with 
the experimental bandstructure, in striking contrast with 
our results. These data further highlight the fact that the 
approximations that have been commonly adopted in the 
literature tend to reduce the systematic, sometimes dra- 
matic, overcorrection of the LDA bandstructure which is 
obtained within a rigorous application of the SIC-LDA 
functional. 



C. Materials with d electrons 

As a first step towards studying the effect of SIC on 
transition metal compounds we begin with the case of an 
isolated d-electron atom in a cubic supercell; this allows 
us to determine the effect of SIC on d states while avoid- 
ing complications arising from bandwidth and ligands. In 
particular, we choose for simplicity the Zn 2+ ion, which 
has a completely filled valence shell. Since the Wannier 
transformation tends to mix wavefunctions that overlap 
in space, it is necessary to include the semicore 3s and 
3p states explicitly as valence orbitals; the Zn 3s, 3p and 
3c? orbitals have important spatial overlap (see Fig. [3]), 
in spite of being far from each other in energy. As a con- 
sequence, the pseudopotentials (Troullier and Martins^, 
with cutoff radius rc— 1 a.u. for all channels) are fairly 
hard and impose a relatively stiff plane wave cutoff of 
180 Ry. We use a cell dimension ao = 16 a.u. which is 
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r (bohr) 


Esic (eV) 


(Vsic) (eV) 


0.5479 


-1.590 


-17.972 


0.5480 


-1.591 


-17.976 


0.6289 


-1.172 


-16.548 



TABLE II: Radius from the nucleus r, SIC energy contribu- 
tion per electron, Esic, and the average of the SIC potential, 
(Vsic) calculated for the three groups of sp d 5 hybrids in the 
isolated Zn 2+ atom. 



FIG. 2: A representative sp 3 d 5 hybrid in the isolated Zn 
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FIG. 3: All-electron and pseudo orbitals for the isolated, neu- 
tral Zn atom, a) 3s; b) 3p; c) 3d. Note the large spatial over- 
lap between s, p and d orbitals, whose maximum is located 
at approximately the same radial distance from the nucleus. 




FIG. 4: Position of the centers of the sp d orbitals. 



large enough so that the spurious crystal-field splitting 
between e g and ti g orbitals is small (5 x 10~ 5 hartree). 

Interestingly, the Wannier localization process of Sec- 
tion [III yields a set of nine similar-looking sp 3 d 5 hy- 
brids (see Fig. [2] for a representative orbital). Upon 
closer inspection of their Wannier centers and their self- 
interaction energies, however, these orbitals differ. In 
fact they are divided into three groups of three members, 
whose main characteristics are summarized in Tab. [TXJ 
It is apparent from the Table that the first two groups 
are practically identical and in fact they form a group 
of six which is artificially broken into two by the tiny 
crystal-field splitting imposed by the cubic symmetry of 
the periodic lattice. The third group, however, is physi- 
cally distinct. To visualize the splitting we indicate the 
centers (i.e. the mean value of the position operator) of 
the orbitals in Fig. dj by highlighting the members of 
the group of six as light blue (darker) small spheres and 
the group of three as white (lighter) small spheres; the 
position of the Zn ion is indicated by a larger sphere. 

Unlike the sp 3 hybrids, which were characterized by 
tctrahcdral symmetry and hence did not cause any split- 
ting in the eigenvalue spectrum within the 3p manifold 
in Ar, the sp 3 d 5 hybrids of Zn are therefore inequivalent. 
This asymmetry is reflected in the eigenvalue spectrum 
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FIG. 5: Schematic diagram of the splittings in the eigenvalue 
spectrum of the SIC-LDA Hamiltonian for the Zn 2+ ion. Val- 
ues are in eV, the spacings are not to scale. 

of the SIC-LDA ground state of the Zn 2+ ion, which is 
represented in Fig.[5j The 3p multiplet, which lies about 
100 eV below the vacuum level, is split into a doublet and 
a singlet by 0.43 eV; however the most dramatic effect is 
found in the 3d electrons, which are extremely important 
for the complex chemistry of the transition metal com- 
pounds. The 3d multiplet is split into three levels (2,2,1) 
by the SIC, with energy separation of 1.12 eV and 0.89 
eV, i.e. a total of 2 eV between the lowest and the highest 
d state. This symmetry breaking and consequent split- 
ting of d levels due to the non-spherical SIC potential 
was already pointed out by Arai and Fujiwara in Ref. 
Il2l . That such a drastically unphysical splitting occurs 
in a spherically symmetric d electron atom renders the 
PZ-SIC formalism unreliable for d electron solids, where 
the interplay of crystal field effects and bandwidth plays 
a dominant role in determining the overall physical prop- 
erties of the compound. 



IV. DISCUSSION 

Our Wannier basis SIC implementation points to two 
problems in the PZ-SIC formalism. The first, the sym- 
metry breaking and splitting of d levels due to the non- 
spherical SIC potential is a problem even for the applica- 
tio of PZ-SIC to atoms. The second, the overestimation 
of the SIE and consequently of splitting in the eigenvalue 
spectrum and band gaps, becomes more acute in many 
electron molecules and solids. Concerning the overesti- 
mation of the gap, we argue that to correct the LDA 
bandstructure, not only the local self-interaction of the 
Wannier charges must be taken into account, but also 



(and especially) the screening properties of the extended 
solid upon electron addition/removal; this physical ingre- 
dient is completely absent in the PZ-SIC-LDA functional, 
which is able to capture the dependency on the environ- 
ment only through the spatial distribution of the Wannier 
charges. This is insufficient for a complete picture: We 
have seen in our examples that the effect of the crystal 
field on the Wannier densities causes remarkably insignif- 
icant variations of the SIC correction to the eigenvalue. 
In particular, for ionic (or rare gas) solids the individual 
constituents are corrected identically to the isolated ion; 
this produces a systematic, gross overestimation of band 
gaps. 

The overestimation of the band gap points to two in- 
teresting and as yet unanswered questions regarding the 
physics of the SIE in many-electron systems: What does 
the self interaction mean in a many-electron system, and 
how does it relate to electronic relaxation? In particu- 
lar, a theory that is self-interaction free in the Koopman 
theorem sense, that is without relaxation corrections, 
will have over-estimated Hartree-Fock band gaps in the 
solid. In fact the self-interaction error is environment- 
dependent, and so the relaxation is not distinct from SIE 
but is intimately related to it. Rigorous theories to in- 
corporate the dependence of the SIE on the dielectric 
screening environment, such as the GW method, tend to 
be costly, even if their range of applicability is steadily 
growing 2 ^. Whenever the problem is an ion embedded 
in a solid with small dispersion and distinct atomic char- 
acter, atomic approximations can be quite effective. For 
example, LDA+U has been used recently as an effective 
technique to cure the SIE, when the value of the Hub- 
bard U parameter is obtained self-consistently within a 
linear-response approach 29 (i.e. it has built-in the dielec- 
tric response of the medium); LDA+U is itself close in 
spirit to the SIC approach, although it was derived from a 
substantially different starting point. However, when the 
covalent character of a given compound is stronger (most 
transition- metal oxides), the reliability of an atomic cor- 
rection applied only to selected bands becomes question- 
able, hence the interest for a more uniform treatment of 
the occupied bands. 

A workaround to this problem within SIC-LDA could 
be to scale down the SIC contribution by a suitable pref- 
actor. This was the approach adopted in the pseudo-SIC 
formalism of Filippetti and Spaldinii where a reduction 
of the atomic SIC by a factor of 0.5 was included to 
account for relaxation effects; within full-SIC, Bylaska, 
Tsemekhman and Gao — found that a factor of 0.2 was 
appropriate to describe defects in Si and C compounds. 
Recent work for molecules^ showed that a scaling fac- 

tor of ( - s — J , where r a is the noninteracting kinetic en- 
ergy density of a spin electrons, and — ^^fey is 
the Weizacker kinetic energy density, gives improved be- 
havior for atoms and molecules. All of these methods 
improve the agreement with the experimental bandgaps, 
while still retaining the main physical advantage of SIC: 
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the Hartree Fock-like treatment of the on-site Coulomb 
interactions. 

Scaling down the SIC, however, docs not remove the 
unphysical symmetry breaking, which is especially seri- 
ous in <i-clcctron (and presumably /-electron) materi- 
als, and is due to the lack of unitary invariance of the 
functional. Therefore we propose that the most promis- 
ing route to incorporating the SIC while preserving uni- 
tary invariance seems to be the use of hybrid functionals, 
which incorporate a fraction of HF exchange. Hybrid 
functionals have yielded very encouraging results for a 
wide class of systems for the bandgap, structural and 
electronic properties; the Wannier function methods pre- 
sented in this work might be useful in devising efficient 
implementations of the Fock exchange within a plane- 
wave pseudopotential formalism. 



V. CONCLUSIONS 

In summary, we have demonstrated that some com- 
mon problems of SIC-LDA, including multiple local min- 
ima and size-consistency issues, can be avoided by lifting 
the spherical approximation. However, our calculations 
expose two other serious drawbacks of SIC-LDA. First, 
we find that the application of the SIC leads to a dra- 
matic overcorrection of the electronic bandgap compared 
to the LDA solution. Second, we point out a worrisome, 
unphysical symmetry breaking of spherically symmetric 
atoms containing d electrons; based on a perturbative 
analysis we argue that this drawback might be a gen- 
eral feature of state-dependent functionals. Our results 
highlight the deficiencies of state-dependent corrections 
to approximate Kohn-Sham theories, and suggest that 
rotationally invariant corrections (such as the Hartree- 
Fock exchange in hybrid functionals) are more promising. 



APPENDIX A: SYMMETRY BREAKING AND 
BOYS LOCALIZATION 

To better understand whether the origin of this break- 
down of spherical symmetry is related to the special fea- 
tures of SIC-LDA or is a more general effect that might 
concern any orbital-dependent functional, it is useful to 
remove all unnecessary complications and look at the ef- 
fect of the simplest possible orbital-dependent perturba- 
tion on the multiplet structure of a spherically symmetric 
atom with a filled d shell. A very practical choice is the 
quadratic spread by Marzari and Vandcrbilt, which is 
better known as Boys quadratic spread in isolated atoms 
and molecules. Working in the framework of perturba- 
tion theory, we start by adding a small contribution to 
the KS total energy: 



First we need to minimize 51 with respect to unitary 
transformations of the occupied orbitals. We start from 
a representation of angular momentum eigenstates corre- 
sponding to the n — 3 shell of Zn, i.e. a total of 9 states 
identified by I and m quantum numbers. It is clear from 
the above equation that the minimum spread is achieved 
by maximizing the second terms in the square bracket 
above (the first is invariant); these are the diagonal ele- 
ments of the 3D position operator. Equivalently we seek 
the transformation that minimizes the off-diagonal ele- 
ments of the three projected position operators, which 
do not commute. The real-space position operators are 
particularly simple to calculate on this basis, e.g. for X 
we have: 

(lm\X\l'm') = J ^ lm {v)x^ m ,{v)d 3 r ; 

x can be written as a real solid harmonic function with 
I = 1: 




so that the above matrix element is simply evaluated in 
terms of a radial integral and a Gaunt coefficient, G: 

{lm\X\l'm')=-sJ— J r^UrWi'm'^drG}*^,. 

Angular momentum selection rules yield a set of three 
sparse matrices with zero diagonal elements (angular mo- 
mentum eigenstates are all centered around the origin) 
that depend on two values only, which are the sp and 
pd radial overlaps; the solution to the problem, apart 
from a trivial scaling factor, is therefore determined by 
a single parameter. We used the values from the LDA 
solution of the isolated Zn pseudoatom to construct the 
three matrices. We then induced a small random unitary 
mixing to break the symmetry, and we further optimized 
the quadratic spread to the minimum until we obtained 
a set of nine sp 3 d 5 localized hybrids. 

By looking at the positions of the centers (radius from 
the origin) and at the individual values of the quadratic 
spread of each hybrid, it is clear that the simplified MV 
spread functional qualitatively reproduces the same lo- 
calization pattern which was induced by SIC-LDA. In 
particular, the orbitals are split into two inequivalcnt 
groups, one of six and one of three members; the much 
smaller splitting of the group of six into two subgroups of 
three (which was observed in our SIC-LDA atomic cal- 
culations) is not reproduced here since in this case we do 
not use the supercell approach. To appreciate the im- 
pact of this perturbation on the Hamiltonian, we now 
take the functional derivative of f2 with respect to the 
wavefunctions; since we are at the variational minimum 
with respect to unitary rotations, this functional deriva- 
tive yields a well-defined Hermitian operator: 



E x = E KS + Xfl, f> = ]T[(r 2 ) 4 -rj] 

i 



H X = Hks + A[i? 2 - 2 ]T fi.(wi\B] 

i 
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The R 2 operator is not orbital dependent, but rather 
a harmonic 3D potential acting on all orbitals which 
preserves spherical symmetry. The second term in the 
square bracket, however, does break the symmetry; in 
particular, it introduces couplings between states be- 
longing to different angular momentum multiplcts. The 
eigenvalues of this operator reproduce qualitatively the 
splittings observed in the SIC-LDA solution, with the 3 
degenerate p states split into 2-1, and the 5 degenerate d 



states split into 2-2-1. Therefore, even in the simple case 
of a harmonic state-dependent perturbation, the spher- 
ical symmetry of the isolated atom is broken, and the 
angular momentum is no longer a good quantum num- 
ber. This is a very undesirable (and artificial) effect that 
poses serious problems for the practical applicability of 
state-dependent potentials to the physics of transition 
metal oxides. 
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